Loading libraries

library(dada2); packageVersion("dada2")
## [1] '1.18.0'
library(ggplot2); packageVersion("ggplot2")
## [1] '3.3.3'
library(phyloseq); packageVersion("phyloseq")
## [1] '1.34.0'
library(data.table); packageVersion("data.table")
## [1] '1.14.0'
library(dplyr)
library(plyr)
library(magrittr)
library(RColorBrewer)
library(picante)
library(igraph)
library(cowplot)
library(plotly)
library(lemon)
library(GGally)
library(randomcoloR)
library(knitr)
library(kableExtra)
library(ggpmisc)
library(ggpubr)
library(decontam)
library(stringr)
library(tibble)

Importing data

SVt = readRDS("../seqtab_final.rds")
# SVt_cnn = collapseNoMismatch(SVt, minOverlap = 20,
#                               identicalOnly = FALSE, vec = TRUE, verbose = TRUE)
# saveRDS(SVt_cnn, "../seqtab_final_cnn.rds")
# SVt = readRDS("../seqtab_final_cnn.rds")
SVt = SVt[!rownames(SVt) %in% "Undetermined", ]

Tax = readRDS("../tax_final.rds")

metadata<-read.csv("../../220319_for_graph.csv", header=TRUE)
# metadata[, 10:47]  <- as.numeric(metadata[, 10:47])
metadata$Snowpack_Layer = gsub("Glacial Ice", "GL ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Sup. Ice", "SUP ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Top", "TOP", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Mid", "MID", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Basal", "BASAL", metadata$Snowpack_Layer)

order<-rownames(SVt)
metadata$Sample_ID=as.character(metadata$Sample_ID)
metadata_ord<-metadata[match(order, metadata$Sample_ID),]
rownames(metadata_ord) <- metadata_ord$Sample_ID
metadata_ord$Sample_Name <- factor(metadata$Sample_Name, levels=unique(metadata$Sample_Name))
ps <- phyloseq(otu_table(SVt, taxa_are_rows=FALSE),
               sample_data(metadata_ord),
               tax_table(Tax))
ps.b <- prune_taxa(taxa_sums(ps) > 0, ps)
ps.b = subset_taxa(ps.b, Kingdom %in% c("Archaea", "Bacteria"))
ps.b = subset_samples(ps.b, Fox_Aspect != "N_NE_D")

Testing batch effects

#remove samples with 0 reads
ps.batch = prune_samples(sample_sums(ps.b)>=1, ps.b)

readsumsdf_samples<-data.frame(nreads = sample_sums(ps.batch),
                               samples = rownames(as.data.frame(sample_sums(ps.batch))),
                               batch = sample_data(ps.batch)$Batch.no.)
readsumsdf_samples_f = readsumsdf_samples[- grep("Control", readsumsdf_samples$batch),]

my_comparisons=list(c("1","2"),c("1","3"),c("1","4"),c("1","5"),
                    c("2","3"), c("2","4"), c("2","5"),
                    c("3","4"), c("3","5"),
                    c("4","5"))

ggplot(readsumsdf_samples_f, aes(x=batch, y=nreads,
                                 color=batch, fill=batch)) +
  geom_violin(alpha=0.3) +
  geom_jitter(shape=16, position = "jitter") +
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Batch Number") +
  guides(fill=FALSE) +
  xlab("Extraction batch") +
  ylab("Number of reads obtained") +
  theme(axis.text.x = element_text(vjust = 1,
                                   hjust = 1,
                                   size = 12),
        strip.text = element_text(size = 12, margin = margin()))

1 - Library sizes

sort(sample_sums(ps.b))
##     C05-cDNA-29     D12-cDNA-MM     C10-cDNA-34 A02-cDNA-MM-H20     C08-cDNA-32 
##               0               0               5               7              25 
##     C07-cDNA-31     B07-cDNA-17 E01-cDNA-MM-H20     B11-cDNA-22     B09-cDNA-19 
##              27              43              55             123             187 
##     D01-cDNA-37     A01-cDNA-MM      A11-cDNA-9     D10-cDNA-46     B10-cDNA-21 
##             262             290            2320            3844            9170 
##      A10-cDNA-8     B05-cDNA-15       B08-DNA-4       B11-DNA-7     B08-cDNA-18 
##           15595           15781           16781           29350           29808 
##     C03-cDNA-27      A08-cDNA-6      A06-cDNA-4      A09-cDNA-7      A07-cDNA-5 
##           31806           33041           37485           40604           44392 
##       A11-DNA-9      A05-cDNA-3       A08-DNA-6     B02-cDNA-12       D08-DNA-8 
##           50662           52562           53872           57398           57493 
##      E02-DNA-MM     C02-cDNA-26     D07-cDNA-43       D06-DNA-6       D02-DNA-2 
##           58537           59570           59805           60653           61230 
##      A01-DNA-MM       B12-DNA-8     B12-cDNA-23      A12-DNA-10       C07-DNA-5 
##           61807           62146           62235           63891           64148 
##     B04-cDNA-14       A07-DNA-5  A02-DNA-MM-H20      A04-cDNA-2     B03-cDNA-13 
##           68029           69142           69499           70613           76895 
##       D09-DNA-9      B02-DNA-12       A05-DNA-3     C04-cDNA-28       C05-DNA-3 
##           80028           83134           84732           86983           87568 
##       C04-DNA-2     D06-cDNA-42       A09-DNA-7     B01-cDNA-11      B04-DNA-14 
##           96798          101335          108915          110464          113222 
##    C03-DNA-B2-1      B03-DNA-13       C08-DNA-6     A12-cDNA-10     B06-cDNA-16 
##          114048          118548          122463          123046          123687 
##     C12-cDNA-36     D02-cDNA-38       D03-DNA-3       C06-DNA-4  E03-DNA-MM-H20 
##          123719          137347          138813          139808          140502 
##       D07-DNA-7      C12-DNA-10       D04-DNA-4       B09-DNA-5     D08-cDNA-44 
##          140668          148596          151009          155211          158902 
##       A06-DNA-4       C01-DNA-9       B10-DNA-6       A10-DNA-8       C09-DNA-7 
##          166098          169126          170463          170498          173577 
##      B01-DNA-11       E01-DNA-4     C06-cDNA-30     C09-cDNA-33     D03-cDNA-39 
##          174897          175874          188016          190708          192008 
##      A03-cDNA-1       A03-DNA-1       B06-DNA-2     C11-cDNA-35       A04-DNA-2 
##          193167          196436          198649          198996          200229 
##       B07-DNA-3     D04-cDNA-40       C10-DNA-8     D11-cDNA-47    D01-DNA-B4-1 
##          204577          204668          212618          220186          238199 
##     D09-cDNA-45       D11-DNA-2       D05-DNA-5    D10-DNA-B5-1     C01-cDNA-25 
##          240324          252560          267031          280585          293592 
##     D05-cDNA-41       D12-DNA-3 
##          302761          312011
df <- as.data.frame(sample_data(ps.b)) # Put sample_data into a ggplot-friendly data.frame
df$LibrarySize <- sample_sums(ps.b)
df <- df[order(df$LibrarySize),]
df$Index <- seq(nrow(df))
ggplot(data=df, aes(x=Index, y=LibrarySize, color=Sample_type)) + geom_point()

2 - Evaluating dataset contamination with decontam

#remove samples with 0 reads
ps.b = prune_samples(sample_sums(ps.b)>=1, ps.b)

sample_data(ps.b)$is.neg <- sample_data(ps.b)$Sample_type == "CONTROL"
contamdf0.5.prev <- isContaminant(ps.b, method="prevalence", threshold = 0.1,
                               neg="is.neg")

ps.final.noncontam <- prune_taxa(!contamdf0.5.prev$contaminant, ps.b)
ps.final.noncontam.noconts = subset_samples(ps.final.noncontam, Sample_type != "CONTROL")

ps.final.noncontam.noconts.nodoubs = subset_samples(ps.final.noncontam.noconts,
    Sample_ID != "D02-DNA-2" & Sample_ID != "B11-DNA-7" & Sample_ID != "C10-cDNA-34" & Sample_ID != "A08-cDNA-6")

ps.final.noncontam.noconts.nodoubs = prune_samples(sample_sums(ps.final.noncontam.noconts.nodoubs)>=1000,
                                                   ps.final.noncontam.noconts.nodoubs)

ps.final <- prune_taxa(taxa_sums(ps.final.noncontam.noconts.nodoubs) > 0, ps.final.noncontam.noconts.nodoubs)
ps.final = prune_samples(sample_sums(ps.final)>=1, ps.final)

ps.final.r <-  transform_sample_counts(ps.final, function(x) x / sum(x))
ps.final.r.f = filter_taxa(ps.final.r, function(x) sum(x) > .001, TRUE)
ps.final.r
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 35534 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 35534 taxa by 6 taxonomic ranks ]
ps.final.r.f
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2886 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 2886 taxa by 6 taxonomic ranks ]
ps.final
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 35534 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 35534 taxa by 6 taxonomic ranks ]
sum(rowSums(otu_table(ps.final)))
## [1] 6750655
gl_ice_spls = as.character(unique(sample_data(ps.final)[str_detect(sample_data(ps.final)$Sample_Name, "Gl ice"),]$Sample_Name))

Setting order for all plots

sample_data(ps.final.r)$Month = factor(sample_data(ps.final.r)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r)$Fox_Aspect = factor(sample_data(ps.final.r)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r)$Snowpack_Layer = factor(sample_data(ps.final.r)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

sample_data(ps.final)$Month = factor(sample_data(ps.final)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final)$Fox_Aspect = factor(sample_data(ps.final)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final)$Snowpack_Layer = factor(sample_data(ps.final)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

sample_data(ps.final.r.f)$Month = factor(sample_data(ps.final.r.f)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r.f)$Fox_Aspect = factor(sample_data(ps.final.r.f)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r.f)$Snowpack_Layer = factor(sample_data(ps.final.r.f)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

3a - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Fox_Aspect, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

ps.final.r.glom <- tax_glom(subset_samples(ps.final.r, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)), taxrank = 'Phylum', NArm = FALSE)
ps.final.r.glom.psdf <- data.table(psmelt(ps.final.r.glom))
ps.final.r.glom.psdf$Phylum <- as.character(ps.final.r.glom.psdf$Phylum)

ps.final.r.glom.psdf[, mean := mean(Abundance, na.rm = TRUE), 
                 by = "Phylum"]
ps.final.r.glom.psdf[(mean <= 0.01), Phylum := "Other"]
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Bacteria"] <- "Unc. Bacteria"
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Archaea"] <- "Unc. Archaea"
#Removing Proteobacteria from previous table
ps.final.r.glom.psdf.noProteo<-ps.final.r.glom.psdf[!grepl("Proteobacteria", ps.final.r.glom.psdf$Phylum),]

#New table, with Proteobacterial classes only
ps.final.r.ProtCl = subset_taxa(ps.final.r, Phylum == "Proteobacteria")
ps.final.r.ProtClglom <- tax_glom(ps.final.r.ProtCl, taxrank = 'Class', NArm = FALSE)
ps.final.r.ProtClglom.psdf <- data.table(psmelt(ps.final.r.ProtClglom))
ps.final.r.ProtClglom.psdf$Class <- as.character(ps.final.r.ProtClglom.psdf$Class)

ps.final.r.ProtClglom.psdf[, mean := mean(Abundance, na.rm = TRUE), 
                       by = "Class"]
ps.final.r.ProtClglom.psdf[(mean <= 0.001), Class := "Other Proteobacteria"]
ps.final.r.ProtClglom.psdf.noP <- subset(ps.final.r.ProtClglom.psdf, select = -Phylum)

names(ps.final.r.glom.psdf.noProteo)[names(ps.final.r.glom.psdf.noProteo) == 'Phylum'] <- 'Taxa'
names(ps.final.r.ProtClglom.psdf.noP)[names(ps.final.r.ProtClglom.psdf.noP) == 'Class'] <- 'Taxa'
# ps.final.r.glom.psdf.noProteo = subset(ps.final.r.glom.psdf.noProteo, select=-c(Class))
ps.final.r.ProteoClAllPhy <- rbind(ps.final.r.glom.psdf.noProteo, ps.final.r.ProtClglom.psdf.noP)

tol18rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455")

ps.final.r.ProteoClAllPhy$Taxa = factor(ps.final.r.ProteoClAllPhy$Taxa, 
              levels = c("Alphaproteobacteria", "Deltaproteobacteria",
                         "Gammaproteobacteria",
                         "Other Proteobacteria",
                         "Actinobacteria",
                         "Acidobacteria","Bacteroidetes",
                         "Chloroflexi","Cyanobacteria",
                         "Firmicutes","Gemmatimonadetes","Other"))
ggplot(ps.final.r.ProteoClAllPhy, 
       aes(x = Month, y = Abundance, fill = Taxa)) + 
  facet_grid(DNA_or_cDNA~Fox_Aspect,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill") +
  scale_fill_manual(values = tol18rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 12),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 12)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance\n")

3b - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Snowpack_Layer

ggplot(ps.final.r.ProteoClAllPhy, 
                                   aes(x = Month, y = Abundance, fill = Taxa)) + 
  facet_grid(DNA_or_cDNA~Snowpack_Layer,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill") +
  scale_fill_manual(values = tol18rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 11),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 12)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance\n")

3c - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Month

tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#AA4455")

ggplot(ps.final.r.ProteoClAllPhy,
                                   aes(x = Sample_Name, y = Abundance, fill = Taxa)) +
  facet_grid(DNA_or_cDNA~Month,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) +
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

4 - Top SVs and their sequences

Top20SVs_names = names(sort(taxa_sums(ps.final.r), TRUE)[1:20])
Top20SVs_names_df=as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6])
rownames(Top20SVs_names_df) = c()
Top20SVs_names_df
##            Phylum               Class                 Order            Family
## 1  Proteobacteria Alphaproteobacteria           Rhizobiales  Beijerinckiaceae
## 2  Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 3  Proteobacteria Alphaproteobacteria           Rhizobiales         Labraceae
## 4  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 5  Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 6  Proteobacteria Alphaproteobacteria           Rhizobiales Xanthobacteraceae
## 7  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 8  Proteobacteria Alphaproteobacteria           Rhizobiales  Beijerinckiaceae
## 9  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 10 Proteobacteria Alphaproteobacteria       Acetobacterales  Acetobacteraceae
## 11 Actinobacteria      Actinobacteria         Micrococcales Microbacteriaceae
## 12 Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 13 Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 14 Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 15 Actinobacteria      Actinobacteria         Micrococcales    Micrococcaceae
## 16 Actinobacteria      Actinobacteria     Corynebacteriales      Nocardiaceae
## 17 Actinobacteria      Actinobacteria            Frankiales              <NA>
## 18 Actinobacteria      Actinobacteria         Micrococcales Cellulomonadaceae
## 19 Actinobacteria      Actinobacteria            Frankiales              <NA>
## 20 Actinobacteria      Actinobacteria            Frankiales              <NA>
##                                         Genus
## 1                                       Bosea
## 2                                 Cupriavidus
## 3                                      Labrys
## 4                                Sphingomonas
## 5  Burkholderia-Caballeronia-Paraburkholderia
## 6                              Bradyrhizobium
## 7                                Sphingomonas
## 8                            Methylobacterium
## 9                                Sphingomonas
## 10                                       <NA>
## 11                            Salinibacterium
## 12                                   Massilia
## 13                               Sphingomonas
## 14                                    Delftia
## 15                               Arthrobacter
## 16                                       <NA>
## 17                                       <NA>
## 18                                       <NA>
## 19                                       <NA>
## 20                                       <NA>
rownames(as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6]))
##  [1] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTGTCCGGGAAGATAATGACTGTACCGGAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGACTCTTAAGTCGGGGGTGAAAGCCCAGGGCTCAACCCTGGAATTGCCTTCGATACTGGGAGTCTTGAGTTCGGAAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGTCCGAAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
##  [2] "TGGGGAATTTTGGACAATGGGGGCAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAATCGCGCTGGTTAATACCTGGCGTGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGACAGGCGTGAAATCCCCGGGCTTAACCTGGGAATTGCGCTTGTGACTGCAAGGCTAGAGTGCGTCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGACGTGACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
##  [3] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGACGGCCTTAGGGTTGTAAAGCTCTTTTAACAGGGACGATAATGACGGTACCTGTAGAATAAGCCCCGGCAAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATTGTTAAGTCGGGGGTGAAATCCTGAGGCTCAACCTCAGAACTGCCTTCGATACTGGCGATCTTGAGTTCGGAAGAGGTTGGTGGAACAGCTAGTGTAGAGGTGAAATTCGTAGATATTAGCTAGAACACCAGTGGCGAAGGCGGCCAACTGGTCCGATACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [4] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCTGGTGCTCAACACCAGAACTGCCTTTTAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [5] "TGGGGAATTTTGGACAATGGGGGAAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAAACCTTTGCACTAATACTGTGAGGGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTCGTTAAGACAGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTTGTGACTGGCGAGCTAGAGTATGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGCCAATACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
##  [6] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTGTGCGGGAAGATAATGACGGTACCGCAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGTGCGTAGGCGGGTCTTTAAGTCAGGGGTGAAATCCTGGAGCTCAACTCCAGAACTGCCTTTGATACTGAAGATCTTGAGTTCGGGAGAGGTGAGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGATACTGACGCTGAGGCACGAAAGCGTGGGGAGCAAACA"                         
##  [7] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [8] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTATCCGGGACGATAATGACGGTACCGGAGGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGCGTTTTAAGTCGGGGGTGAAAGCCTGTGGCTCAACCACAGAATGGCCTTCGATACTGGGACGCTTGAGTATGGTAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGACCATCACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
##  [9] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTACCCGGGAAGATAATGACTGTACCGGGAGAATAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTCAGAGGTGAAAGCCTGGAGCTCAACTCCAGAACTGCCTTTGAGACTGCATCGCTTGAATCCAGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
## [10] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGTCTTCGGATTGTAAAGCACTTTCGGCGGGGACGATGATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATGACTGGGCGTAAAGGGCGCGTAGGCGGTTGCATAAGTTAGATGTGAAATTCCCGGGCTTAACCTGGGGACTGCATTTGATACTATGTGGCTTGAGTATGGAAGAGGGTCGTGGAATTCCCAGTGTAGAGGTGAAATTCGTAGATATTGGGAAGAACACCGGTGGCGAAGGCGGCGACCTGGTCCATAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
## [11] "TGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCAACGCCGCGTGAGGGACGACGGCCTTCGGGTTGTAAACCTCTTTTAGTAGGGAAGAAGCGAAAGTGACGGTACCTGCAGAAAAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTGGGGGCTCAACCCCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCAATGGCGAAGGCAGATCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"                    
## [12] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGAAGAAGGCCTTCGGGTTGTAAAGCTCTTTTGTCAGGGAAGAAACGGTGTGGGCTAATATCCTGCACTAATGACGGTACCTGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGTCTGTCGTGAAATCCCCGGGCTCAACCTGGGAATTGCGATGGAGACTGCAAGGCTAGAATCTGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGTCAAGATTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [13] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTTAGAGGTGAAAGCCTGGAGCTCAACTCCAGAATTGCCTTTAAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
## [14] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGCAGGATGAAGGCCTTCGGGTTGTAAACTGCTTTTGTACGGAACGAAAAAGCTCCTTCTAATACAGGGGGCCCATGACGGTACCGTAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTATGTAAGACAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTTGTGACTGCATGGCTAGAGTACGGTAGAGGGGGATGGAATTCCGCGTGTAGCAGTGAAATGCGTAGATATGCGGAGGAACACCGATGGCGAAGGCAATCCCCTGGACCTGTACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [15] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGTAGGGAACAAGGCCACACGTGGTGTGGTTGAGGGTACTTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCCGTGAAAGTCCGGGGCTCAACTCCGGATCTGCGGTGGGTACGGGCAGACTAGAGTGATGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCATTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"         
## [16] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAAGCGAGAGTGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTGTCCGGAATTACTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCACGTCGGCTGTGAAAACCCATCGCTCAACGGTGGGCTTGCAGTCGATACGGGCTGACTTGAGTACTGCAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"                    
## [17] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"                        
## [18] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGGAGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGAAGAAGTGGTCGGGTTCTCTCGGTCGTGACGGTACCTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTCGAGGCTCAACCTCGGGCTTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCCGCAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"      
## [19] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCAAACA"                        
## [20] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAACACAGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTCTGTCGCGTCGGCTGTGAAAACTCGGGGCTCAACCCCGAGCCTGCAGTCGATACGGGCAGACTAGAGTGCTGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"

5a - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (no filter)

5b - Genus-level stacked barchart by DNA_or_cDNA~Snowpack_Layer (filter)

Oh good. But this need improving. Getting the unclassified component of the families for the top main taxa in it, but without a guarantee for them popping into it. Also, basing it on the >0.1% dataset as per Arwyn’s advice.

ps.final.r.glom.genus.psdf_2 <- data.table(speedyseq::psmelt(ps.final.r.glom.genus))
ps.final.r.glom.genus.psdf_2$Genus <- as.character(ps.final.r.glom.genus.psdf_2$Genus)

ps.final.r.glom.genus.psdf_2[, mean := mean(Abundance, na.rm = TRUE), 
                 by = "Genus"]

ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Sphingomonadaceae"] <- "Unc. Sphingomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Burkholderiaceae"] <- "Unc. Burkholderiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Labraceae"] <- "Unc. Labraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Acetobacteraceae"] <- "Unc. Acetobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Microbacteriaceae"] <- "Unc. Microbacteriaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Deinococcaceae"] <- "Unc. Deinococcaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Gemmatimonaceae"] <- "Unc. Gemmatimonaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Intrasporangiaceae"] <- "Unc. Intrasporangiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Nocardiaceae"] <- "Unc. Nocardiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Ktedonobacteraceae"] <- "Unc. Ktedonobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Cellulomonadaceae"] <- "Unc. Cellulomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Beijerinckiaceae"] <- "Unc. Beijerinckiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          is.na(ps.final.r.glom.genus.psdf_2$Family)] <- "Unclassified"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &  is.na(ps.final.r.glom.genus.psdf_2$Family) & ps.final.r.glom.genus.psdf_2$Order=="Frankiales"] <- "Unc. Frankiales"

ps.final.r.glom.genus.psdf_2[(mean <= 0.01), Genus := "Other (<1%)"]

ps.final.r.glom.genus.psdf_2$Genus = factor(ps.final.r.glom.genus.psdf_2$Genus,
                         levels = c("Acidiphilium","Arthrobacter",
                                    "Bacillus",
                                    "Bosea",
                                  "Burkholderia-Caballeronia-Paraburkholderia",
                                    "Cryobacterium","Cupriavidus",
                                    "Delftia","Hymenobacter",
                                    "Labrys","Leptothrix",
                                    "Massilia","Methylobacterium",
                                    "Noviherbaspirillum","Oryzihumus",
                                    "Polymorphobacter",
                                    "Salinibacterium",
                                    "Sphingomonas","Other (<1%)"))
tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477",
               "#4477AA", "#77AADD", "#117777", "#44AAAA", 
               "#77CCCC", "#777711", "#AAAA44", "#DDDD77", 
               "#774411", "#AA7744", "#DDAA77", "#771122", 
               "#AA4455", "#DD7788", "azure3")

ggplot(ps.final.r.glom.genus.psdf_2, 
       aes(x = Month, y = Abundance, fill = Genus)) + 
  facet_grid(DNA_or_cDNA~Snowpack_Layer,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

5c - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (filter)

ggplot(ps.final.r.glom.genus.psdf_2, 
                                   aes(x = Month, y = Abundance, fill = Genus)) + 
  facet_grid(DNA_or_cDNA~Fox_Aspect,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

Lest us not forget that size of Other (<1%) corresponds to the different samples included in that category.

I don’t like this, but plotting all samples now, without grouping.

# ps.final.r.glom.genus.psdf_2$Sample_Name <- as.character(ps.final.r.glom.genus.psdf_2$Sample_Name)
# ps.final.r.glom.genus.psdf_2$Sample_Name <- factor(ps.final.r.glom.genus.psdf_2$Sample_Name,
#                                                  levels=unique(ps.final.r.glom.genus.psdf_2$Sample_Name))
ggplot(ps.final.r.glom.genus.psdf_2,
                                   aes(x = Sample_Name, y = Abundance, fill = Genus)) +
  facet_grid(DNA_or_cDNA~Month,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) +
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

6 - Unconstrained ordination

nmds_shapes = c(15,16,17,18,6)
nmdscolors_by_month = c("April" = "#A6048B",
                        "June" = "#0704A6",
                        "Early July" = "#727A0B",
                        "Late July" = "#6C0504")
nmdscolors_by_snpl = c("Top" = "#151AC6",
                       "Mid" = "#C6156D",
                       "Basal" = "#890909",
                       "Sup. Ice" = "#091D89",
                       "Glacial Ice" = "#580303")

ps.final.r.dna = subset_samples(ps.final.r, DNA_or_cDNA != "cDNA")
ps.final.r.cdna = subset_samples(ps.final.r, DNA_or_cDNA != "DNA")

ps.final.r.dna.nmds  <- ordinate(ps.final.r.dna,
                         "NMDS",
                         distance="jsd",
                         maxit = 1e3)
## Run 0 stress 0.2107623 
## Run 1 stress 0.214947 
## Run 2 stress 0.2149514 
## Run 3 stress 0.2243087 
## Run 4 stress 0.2176892 
## Run 5 stress 0.2118836 
## Run 6 stress 0.2034066 
## ... New best solution
## ... Procrustes: rmse 0.1219555  max resid 0.4302237 
## Run 7 stress 0.2150477 
## Run 8 stress 0.2102508 
## Run 9 stress 0.2091659 
## Run 10 stress 0.2107622 
## Run 11 stress 0.2171913 
## Run 12 stress 0.2164589 
## Run 13 stress 0.2098655 
## Run 14 stress 0.2034062 
## ... New best solution
## ... Procrustes: rmse 0.0002446002  max resid 0.001304558 
## ... Similar to previous best
## Run 15 stress 0.2126387 
## Run 16 stress 0.228024 
## Run 17 stress 0.2254012 
## Run 18 stress 0.212504 
## Run 19 stress 0.2184797 
## Run 20 stress 0.2179839 
## *** Solution reached
ps.final.r.cdna.nmds  <- ordinate(ps.final.r.cdna,
                         "NMDS",
                         distance="jsd",
                         maxit = 1e3)
## Run 0 stress 0.1654426 
## Run 1 stress 0.1628229 
## ... New best solution
## ... Procrustes: rmse 0.04769025  max resid 0.162577 
## Run 2 stress 0.1641699 
## Run 3 stress 0.1653273 
## Run 4 stress 0.1724367 
## Run 5 stress 0.1659068 
## Run 6 stress 0.1640257 
## Run 7 stress 0.1649862 
## Run 8 stress 0.1604191 
## ... New best solution
## ... Procrustes: rmse 0.0849541  max resid 0.3370751 
## Run 9 stress 0.1713025 
## Run 10 stress 0.166409 
## Run 11 stress 0.1681634 
## Run 12 stress 0.1650155 
## Run 13 stress 0.1754389 
## Run 14 stress 0.1598542 
## ... New best solution
## ... Procrustes: rmse 0.06484638  max resid 0.2497052 
## Run 15 stress 0.1740954 
## Run 16 stress 0.1622533 
## Run 17 stress 0.1643391 
## Run 18 stress 0.1653448 
## Run 19 stress 0.1609927 
## Run 20 stress 0.1654456 
## *** No convergence -- monoMDS stopping criteria:
##     19: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin
ps.final.r.dna.nmds.plot <- plot_ordination(ps.final.r.dna, ps.final.r.dna.nmds) +
               geom_point(aes(colour=Month, shape=Snowpack_Layer),
                          size=5, alpha = 0.5, stroke=1.5) +
               scale_shape_manual(values=nmds_shapes) +
               scale_colour_manual(values=nmdscolors_by_month) +
               scale_fill_manual(values=nmdscolors_by_month) +
               stat_ellipse(aes(fill=Month),
                          geom = "polygon",
                          type="t",
                          alpha=0.15,
                          linetype = 2,
                          show.legend = FALSE) +
  ggtitle("DNA nMDS") +
  guides(colour = guide_legend(title="Month"))

ps.final.r.cdna.nmds.plot <- plot_ordination(ps.final.r.cdna, ps.final.r.cdna.nmds) +
               geom_point(aes(colour=Month, shape=Snowpack_Layer),
                          size=5, alpha = 0.7, stroke=1.5) +
               scale_shape_manual(values=nmds_shapes) +
               scale_colour_manual(values=nmdscolors_by_month) +
               scale_fill_manual(values=nmdscolors_by_month) +
               stat_ellipse(aes(fill=Month),
                          geom = "polygon",
                          type="t",
                          alpha=0.2,
                          linetype = 2,
                          show.legend = FALSE)+
  ggtitle("cDNA nMDS") +
  guides(colour = guide_legend(title="Month"))
plot_grid(ps.final.r.dna.nmds.plot + theme_bw(), 
          ps.final.r.cdna.nmds.plot + theme_bw())

7 - Alpha-diversity metrics.

month_comparisons = list( c("April", "June"), c("April", "Early July"), 
                          c("April", "Late July"),
                          c("June", "Early July"), c("June", "Late July"),
                          c("Early July", "Late July"))

site_comparisons = list( c("NW_AWS", "N_NE"),
                         c("NW_AWS", "SWS1SE"),
                         c("N_NE", "SWS1SE"))

slayer_comparisons = list( c("TOP", "MID"), c("TOP", "BASAL"), c("TOP", "SUP ICE"), 
                           c("TOP", "GL ICE"),
                           c("MID", "BASAL"), c("MID", "SUP ICE"), 
                           c("MID", "GL ICE"),
                           c("BASAL", "SUP ICE"), c("BASAL", "GL ICE"),
                           c("SUP ICE", "GL ICE"))

All measures for Sample_ID, no stats, facet by DNA/cDNA

plot_richness(ps.final,
              x="Sample_ID", color="Sample_ID",
              measures=c("Observed","Shannon","Simpson")) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Sample_ID") +
  guides(color=guide_legend(ncol=2)) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures by Sample_ID

mini_metadata = metadata %>% 
  select(c(Sample_ID, Sample_Name, Month, DNA_or_cDNA, Fox_Aspect, Snowpack_Layer))

ps.final.measures = estimate_richness(ps.final, measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID))

full_join(ps.final.measures, mini_metadata, 
                                       by="Sample_ID") %>% 
  filter(!is.na(Observed)) %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = T, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Sample_ID Observed Shannon Simpson Sample_Name Month DNA_or_cDNA Fox_Aspect Snowpack_Layer
A03-cDNA-1 286 0.7114 0.2103 April_N_NE_Basal_R April cDNA N_NE BASAL
A03-DNA-1 138 3.2774 0.9407 June_N_NE_Top20_D June DNA N_NE TOP
A04-cDNA-2 163 2.0361 0.7852 April_NW_AWS_Top_R April cDNA NW_AWS TOP
A04-DNA-2 75 1.6409 0.6083 April_ NW_AWS_ Basal_D April DNA NW_AWS BASAL
A05-cDNA-3 515 2.2067 0.7048 June_NW_AWS_Basal_R June cDNA NW_AWS BASAL
A05-DNA-3 1390 5.1230 0.9813 30 July_ N_NE_ Gl ice_D Late July DNA N_NE GL ICE
A06-cDNA-4 125 1.3510 0.6058 June_N_NE_Basal_R June cDNA N_NE BASAL
A06-DNA-4 210 3.9337 0.9650 June_ SWS1SE_ Top20_D June DNA SWS1SE TOP
A07-cDNA-5 1635 4.6591 0.9488 9 July_NW_AWS_Mid_R Early July cDNA NW_AWS MID
A07-DNA-5 609 3.3602 0.8726 9 July_ NW_AWS_ Sup ice_D Early July DNA NW_AWS SUP ICE
A08-DNA-6 563 4.0141 0.8786 April_ NW_AWS_ Mid_D April DNA NW_AWS MID
A09-cDNA-7 1875 5.5902 0.9807 30 July_NW_AWS_Top20_R Late July cDNA NW_AWS TOP
A10-cDNA-8 842 4.7339 0.9520 30 July_N_NE_Gl ice_R Late July cDNA N_NE GL ICE
A10-DNA-8 179 3.3916 0.9493 April_NW_AWS_Top20_D April DNA NW_AWS TOP
A11-cDNA-9 3 0.2518 0.1121 April_NW_AWS_Mid_R April cDNA NW_AWS MID
A11-DNA-9 1819 5.3121 0.9764 9 July_SWS1SE_Gl ice 1_D Early July DNA SWS1SE GL ICE
A12-cDNA-10 54 1.4929 0.6133 June_NW_AWS_Top20_R June cDNA NW_AWS TOP
A12-DNA-10 2719 5.7179 0.9837 30 July_N_NE_Top20_D Late July DNA N_NE TOP
B01-cDNA-11 453 4.0612 0.9622 30 July_NW_AWS_Sup ice_R Late July cDNA NW_AWS SUP ICE
B01-DNA-11 187 3.5277 0.9546 June_NW_AWS_Mid_D June DNA NW_AWS MID
B02-cDNA-12 733 4.5851 0.9683 30 July_N_NE_Top20_R Late July cDNA N_NE TOP
B02-DNA-12 310 2.1151 0.7154 9 July_N_NE_Mid_D Early July DNA N_NE MID
B03-cDNA-13 38 1.1241 0.4091 June_SWS1SE_Basal_R June cDNA SWS1SE BASAL
B03-DNA-13 457 4.2859 0.9559 April_SWS1SE_Mid_D April DNA SWS1SE MID
B04-cDNA-14 67 1.1867 0.4139 9 July_NW_AWS_Top20_R Early July cDNA NW_AWS TOP
B05-cDNA-15 1307 5.7881 0.9881 9 July_NW_AWS_Gl ice 1_R Early July cDNA NW_AWS GL ICE
B06-cDNA-16 134 1.4982 0.6595 April_SWS1SE_Top20_R April cDNA SWS1SE TOP
B06-DNA-2 156 0.6586 0.2271 April_SWS1SE_Basal_D April DNA SWS1SE BASAL
B07-DNA-3 1302 5.3971 0.9910 June_NW AWS_Top20_D June DNA NW_AWS TOP
B08-cDNA-18 21 0.9250 0.5154 June_NW_AWS_Mid_R June cDNA NW_AWS MID
B08-DNA-4 986 5.2137 0.9805 30 July_SWS1SE_Gl ice_D Late July DNA SWS1SE GL ICE
B09-DNA-5 105 2.9238 0.9150 June_N_NE_Mid_D June DNA N_NE MID
B10-cDNA-21 400 4.5700 0.9771 9 July_SWS1SE_Sup ice_R Early July cDNA SWS1SE SUP ICE
B10-DNA-6 377 4.1547 0.9598 June_SWS1SE_Basal_D June DNA SWS1SE BASAL
B12-cDNA-23 29 0.9625 0.3511 30 July_SWS1SE_Slush_R Late July cDNA SWS1SE TOP
B12-DNA-8 1998 5.3382 0.9817 9 July_ NW_AWS_Gl ice 1_D Early July DNA NW_AWS GL ICE
C01-DNA-9 165 3.3370 0.9306 June_SWS1SE_Mid_D June DNA SWS1SE MID
C02-cDNA-26 907 4.1753 0.9419 30 July_SWS1SE_Gl ice_R Late July cDNA SWS1SE GL ICE
C03-cDNA-27 537 3.3393 0.8840 9 July_N_NE_Sup ice_R Early July cDNA N_NE SUP ICE
C03-DNA-B2-1 677 5.1429 0.9887 June_NW_AWS_Basal_D June DNA NW_AWS BASAL
C04-cDNA-28 307 3.3747 0.8749 30 July_NW_AWS_Gl ice_R Late July cDNA NW_AWS GL ICE
C04-DNA-2 715 4.2502 0.9430 30 July_NW_AWS_Gl ice_D Late July DNA NW_AWS GL ICE
C05-DNA-3 1739 5.2347 0.9838 30 July_NW_AWS_Sup ice_D Late July DNA NW_AWS SUP ICE
C06-DNA-4 236 3.4402 0.9006 April_SWS1SE_Top20_D April DNA SWS1SE TOP
C07-DNA-5 2789 5.9749 0.9894 9 July_NW_AWS_Mid_D Early July DNA NW_AWS MID
C08-DNA-6 307 4.5001 0.9837 June_N_NE_Basal_D June DNA N_NE BASAL
C09-cDNA-33 36 1.3284 0.5802 June_N_NE_Top20_R June cDNA N_NE TOP
C09-DNA-7 148 3.8243 0.9690 9 July_N_NE_Gl ice_D Early July DNA N_NE GL ICE
C10-DNA-8 383 3.2676 0.8801 9 July_ SWS1SE_Mid_D Early July DNA SWS1SE MID
C11-cDNA-35 65 1.5747 0.6657 9 July_SWS1SE_Mid_R Early July cDNA SWS1SE MID
C12-cDNA-36 19 0.4213 0.1462 9 July_N_NE_Gl ice_R Early July cDNA N_NE GL ICE
C12-DNA-10 396 2.2692 0.7492 9 July_N_NE_Sup ice_D Early July DNA N_NE SUP ICE
D01-DNA-B4-1 147 0.6873 0.2422 30 July_SWS1SE_Slush_D Late July DNA SWS1SE TOP
D02-cDNA-38 128 2.1030 0.7780 April_SWS1SE_Mid_R April cDNA SWS1SE MID
D03-cDNA-39 266 3.4007 0.9035 April_N_NE_Mid_R April cDNA N_NE MID
D04-cDNA-40 79 1.5083 0.6114 June_SWS1SE_Top20_R June cDNA SWS1SE TOP
D05-cDNA-41 154 1.8281 0.7377 June_N_NE_Mid_R June cDNA N_NE MID
D05-DNA-5 90 0.6634 0.2169 9 July_N_NE_Top20_D Early July DNA N_NE TOP
D06-cDNA-42 841 3.6042 0.9024 9 July_NW_AWS_Sup ice_R Early July cDNA NW_AWS SUP ICE
D06-DNA-6 2292 5.5649 0.9810 9 July_SWS1SE_Sup ice_D Early July DNA SWS1SE SUP ICE
D07-cDNA-43 1444 4.5026 0.9489 9 July_SWS1SE_Gl ice 2_R Early July cDNA SWS1SE GL ICE
D07-DNA-7 467 1.8746 0.5308 9 July_NW_AWS_Top20_D Early July DNA NW_AWS TOP
D08-cDNA-44 120 2.6011 0.8428 9 July_N_NE_Mid_R Early July cDNA N_NE MID
D08-DNA-8 1045 4.1510 0.9245 30 July_NW_AWS_Top20_D Late July DNA NW_AWS TOP
D09-DNA-9 2986 5.9104 0.9895 9 July_SWS1SE_Top20_D Early July DNA SWS1SE TOP

All measures for Month, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Month, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

plot_richness(subset_samples(ps.final, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

Table with means

estimate_richness(subset_samples(ps.final, 
     Snowpack_Layer != "GL ICE" & 
     Sample_Name != "30 July_SWS1SE_Slush_R" &
     Sample_Name != "30 July_SWS1SE_Slush_D" &
     !(Sample_Name %in% gl_ice_spls)), 
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Fox_Aspect, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson Month DNA_or_cDNA
169.2000 1.7006 0.5578 April cDNA
270.7500 3.8497 0.9548 June DNA
277.6667 2.9052 0.7533 April DNA
158.6000 1.5037 0.5823 June cDNA
695.6000 3.5495 0.8756 Early July cDNA
1279.0000 3.7921 0.8385 Early July DNA
785.6667 3.5380 0.7647 Late July cDNA
2229.0000 5.4763 0.9837 Late July DNA

All measures for Month, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)

mini_month_comparisons = list(c("April", "June"), 
                          c("April", "June"))

plot_richness(subset_samples(ps.final, Month != "Early July" & Month != "Late July"),
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

Table with means

estimate_richness(subset_samples(ps.final, 
                                 Month != "Early July" & Month != "Late July"),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Fox_Aspect, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson Month DNA_or_cDNA
163.3333 1.6669 0.5748 April cDNA
385.3333 4.0216 0.9588 June DNA
277.6667 2.9052 0.7533 April DNA
127.7500 1.4706 0.5972 June cDNA

All measures for Fox_Aspect, no stats, facet by DNA/cDNA

#now without stats, only shapes
plot_richness(ps.final,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha=0.5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

Table with means

estimate_richness(ps.final,
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson Month DNA_or_cDNA Fox_Aspect
276.0000 2.0560 0.5569 April cDNA N_NE
183.3333 3.5671 0.9465 June DNA N_NE
83.0000 1.1439 0.4487 April cDNA NW_AWS
272.3333 3.0155 0.8120 April DNA NW_AWS
196.6667 1.5415 0.6112 June cDNA NW_AWS
2054.5000 5.4205 0.9825 Late July DNA N_NE
105.0000 1.5025 0.6412 June cDNA N_NE
250.6667 3.8084 0.9518 June DNA SWS1SE
962.5000 3.8095 0.8133 Early July cDNA NW_AWS
1465.7500 4.1370 0.8436 Early July DNA NW_AWS
878.3333 4.3420 0.9393 Late July cDNA NW_AWS
787.5000 4.6595 0.9602 Late July cDNA N_NE
1870.0000 5.0137 0.9568 Early July DNA SWS1SE
722.0000 4.6893 0.9781 June DNA NW_AWS
236.0000 2.2180 0.6626 Early July DNA N_NE
58.5000 1.3162 0.5102 June cDNA SWS1SE
283.0000 2.7949 0.6945 April DNA SWS1SE
131.0000 1.8006 0.7188 April cDNA SWS1SE
566.5000 2.9505 0.6113 Late July DNA SWS1SE
636.3333 3.5491 0.8639 Early July cDNA SWS1SE
468.0000 2.5689 0.6465 Late July cDNA SWS1SE
225.3333 2.1206 0.6244 Early July cDNA N_NE
1166.3333 4.5453 0.9504 Late July DNA NW_AWS

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

plot_richness(subset_samples(ps.final, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha=0.5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

Table with means

estimate_richness(subset_samples(ps.final, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson Month DNA_or_cDNA Fox_Aspect
276.0000 2.0560 0.5569 April cDNA N_NE
183.3333 3.5671 0.9465 June DNA N_NE
83.0000 1.1439 0.4487 April cDNA NW_AWS
272.3333 3.0155 0.8120 April DNA NW_AWS
515.0000 2.2067 0.7048 June cDNA NW_AWS
80.5000 1.3397 0.5930 June cDNA N_NE
250.6667 3.8084 0.9518 June DNA SWS1SE
1238.0000 4.1317 0.9256 Early July cDNA NW_AWS
1288.3333 3.7366 0.7976 Early July DNA NW_AWS
1164.0000 4.8257 0.9715 Late July cDNA NW_AWS
2719.0000 5.7179 0.9837 Late July DNA N_NE
432.0000 4.3353 0.9716 June DNA NW_AWS
353.0000 2.1921 0.7323 Early July DNA N_NE
58.5000 1.3162 0.5102 June cDNA SWS1SE
283.0000 2.7949 0.6945 April DNA SWS1SE
232.5000 3.0723 0.8214 Early July cDNA SWS1SE
29.0000 0.9625 0.3511 Late July cDNA SWS1SE
537.0000 3.3393 0.8840 Early July cDNA N_NE
1739.0000 5.2347 0.9838 Late July DNA NW_AWS
1887.0000 4.9143 0.9502 Early July DNA SWS1SE
128.0000 2.1030 0.7780 April cDNA SWS1SE

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
                             Month != "Early July" & Month != "Late July"),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means

estimate_richness(subset_samples(ps.final, 
                                 Snowpack_Layer != "GL ICE" &
                             Month != "Early July" & Month != "Late July"),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson Month DNA_or_cDNA Fox_Aspect
276.0000 2.0560 0.5569 April cDNA N_NE
183.3333 3.5671 0.9465 June DNA N_NE
83.0000 1.1439 0.4487 April cDNA NW_AWS
272.3333 3.0155 0.8120 April DNA NW_AWS
196.6667 1.5415 0.6112 June cDNA NW_AWS
105.0000 1.5025 0.6412 June cDNA N_NE
250.6667 3.8084 0.9518 June DNA SWS1SE
722.0000 4.6893 0.9781 June DNA NW_AWS
58.5000 1.3162 0.5102 June cDNA SWS1SE
283.0000 2.7949 0.6945 April DNA SWS1SE
131.0000 1.8006 0.7188 April cDNA SWS1SE

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July")),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means

estimate_richness(subset_samples(ps.final, 
                                 Month == c("Early July", "Late July")),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
Observed Shannon Simpson Month DNA_or_cDNA Fox_Aspect
2054.5000 5.4205 0.9825 Late July DNA N_NE
847.6667 3.1500 0.7550 Early July cDNA NW_AWS
1875.0000 5.5902 0.9807 Late July cDNA NW_AWS
922.0000 4.5363 0.9630 Early July cDNA SWS1SE
907.0000 4.1753 0.9419 Late July cDNA SWS1SE
225.3333 2.1206 0.6244 Early July cDNA N_NE
880.0000 4.2006 0.9338 Late July DNA NW_AWS
2789.0000 5.9749 0.9894 Early July DNA NW_AWS
1684.5000 4.5890 0.9348 Early July DNA SWS1SE

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed."

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July"), 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means

estimate_richness(subset_samples(ps.final, Month == c("Early July", "Late July"), 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
Observed Shannon Simpson Month DNA_or_cDNA Fox_Aspect
2054.5000 5.4205 0.9825 Late July DNA N_NE
847.6667 3.1500 0.7550 Early July cDNA NW_AWS
1875.0000 5.5902 0.9807 Late July cDNA NW_AWS
922.0000 4.5363 0.9630 Early July cDNA SWS1SE
907.0000 4.1753 0.9419 Late July cDNA SWS1SE
225.3333 2.1206 0.6244 Early July cDNA N_NE
880.0000 4.2006 0.9338 Late July DNA NW_AWS
2789.0000 5.9749 0.9894 Early July DNA NW_AWS
1684.5000 4.5890 0.9348 Early July DNA SWS1SE

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Early July” only, “GL ICE” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
                               Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means

estimate_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
                               Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson Month DNA_or_cDNA Fox_Aspect
1238.000 4.1317 0.9256 Early July cDNA NW_AWS
1288.333 3.7366 0.7976 Early July DNA NW_AWS
353.000 2.1921 0.7323 Early July DNA N_NE
232.500 3.0723 0.8214 Early July cDNA SWS1SE
537.000 3.3393 0.8840 Early July cDNA N_NE
1887.000 4.9143 0.9502 Early July DNA SWS1SE

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Late July” only, “GL ICE”, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & Month == "Late July" &
                               Sample_Name != "30 July_SWS1SE_Slush_R" &
                               Sample_Name != "30 July_SWS1SE_Slush_D" & 
                               !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

Table with means

estimate_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & Month == "Late July" &
                               Sample_Name != "30 July_SWS1SE_Slush_R" &
                               Sample_Name != "30 July_SWS1SE_Slush_D" & 
                               !(Sample_Name %in% gl_ice_spls)),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Snowpack_Layer)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(Month, DNA_or_cDNA, Fox_Aspect) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson Month DNA_or_cDNA Fox_Aspect
1164 4.8257 0.9715 Late July cDNA NW_AWS
2719 5.7179 0.9837 Late July DNA N_NE
29 0.9625 0.3511 Late July cDNA SWS1SE
1739 5.2347 0.9838 Late July DNA NW_AWS

All measures for Snowpack_Layer, no stats, facet by DNA/cDNA

#now without stats, only shapes
plot_richness(ps.final,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Snowpack_Layer, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

Table with means

estimate_richness(ps.final,
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(DNA_or_cDNA, Snowpack_Layer) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson DNA_or_cDNA Snowpack_Layer
241.0000 1.3483 0.4825 cDNA BASAL
865.3636 3.4950 0.7849 DNA TOP
352.2222 2.2432 0.6626 cDNA TOP
318.4000 3.2194 0.7535 DNA BASAL
1176.0000 4.8436 0.9720 DNA GL ICE
299.0000 2.1679 0.6880 cDNA MID
1259.0000 4.1072 0.8966 DNA SUP ICE
619.8750 3.6808 0.9025 DNA MID
804.3333 3.8327 0.8087 cDNA GL ICE
557.7500 3.8937 0.9314 cDNA SUP ICE

All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA

plot_richness(subset_samples(ps.final, Month = c("Early July", "Late July")),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

Table with means

estimate_richness(subset_samples(ps.final, 
                                 Month = c("Early July", "Late July")),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(DNA_or_cDNA, Snowpack_Layer) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson DNA_or_cDNA Snowpack_Layer
241.0000 1.3483 0.4825 cDNA BASAL
865.3636 3.4950 0.7849 DNA TOP
352.2222 2.2432 0.6626 cDNA TOP
318.4000 3.2194 0.7535 DNA BASAL
1176.0000 4.8436 0.9720 DNA GL ICE
299.0000 2.1679 0.6880 cDNA MID
1259.0000 4.1072 0.8966 DNA SUP ICE
619.8750 3.6808 0.9025 DNA MID
804.3333 3.8327 0.8087 cDNA GL ICE
557.7500 3.8937 0.9314 cDNA SUP ICE

All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July") &
                                       Snowpack_Layer != "GL ICE" & 
                                       !(Sample_Name %in% gl_ice_spls) &
                                       Sample_Name != c("30 July_SWS1SE_Slush_R",
                                                        "30 July_SWS1SE_Slush_D")),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning in `!=.default`(Sample_Name, c("30 July_SWS1SE_Slush_R", "30
## July_SWS1SE_Slush_D")): longer object length is not a multiple of shorter object
## length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

Table with means

estimate_richness(subset_samples(ps.final, Month == c("Early July", "Late July") &
                                       Snowpack_Layer != "GL ICE" & 
                                       !(Sample_Name %in% gl_ice_spls) &
                                       Sample_Name != c("30 July_SWS1SE_Slush_R",
                                                        "30 July_SWS1SE_Slush_D")),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(DNA_or_cDNA, Snowpack_Layer) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning in `!=.default`(Sample_Name, c("30 July_SWS1SE_Slush_R", "30
## July_SWS1SE_Slush_D")): longer object length is not a multiple of shorter object
## length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
Observed Shannon Simpson DNA_or_cDNA Snowpack_Layer
1635.0000 4.6591 0.9488 cDNA MID
1875.0000 5.5902 0.9807 cDNA TOP
2250.0000 5.2598 0.9659 DNA TOP
592.6667 3.8378 0.9212 cDNA SUP ICE
1586.0000 4.6213 0.9348 DNA MID

All measures for Snowpack_Layer in “Early July”, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & 
                               Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

Table with means

estimate_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & 
                               Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(DNA_or_cDNA, Snowpack_Layer) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson DNA_or_cDNA Snowpack_Layer
850.0000 3.1169 0.8073 cDNA MID
1099.0000 3.7314 0.8676 DNA SUP ICE
1160.6667 3.7859 0.8616 DNA MID
592.6667 3.8378 0.9212 cDNA SUP ICE
1726.5000 3.8925 0.7601 DNA TOP

All measures for Snowpack_Layer in Late July, stats included, facet by DNA/cDNA, “GL ICE” removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_R” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & 
                               Sample_Name != "30 July_SWS1SE_Slush_R" &
                               Sample_Name != "30 July_SWS1SE_Slush_D" &
                               Month == "Late July" & !(Sample_Name %in% gl_ice_spls)),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

Table with means

estimate_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & 
                               Sample_Name != "30 July_SWS1SE_Slush_R" &
                               Sample_Name != "30 July_SWS1SE_Slush_D" &
                               Month == "Late July" & !(Sample_Name %in% gl_ice_spls)),
        measures = c("Observed","Shannon","Simpson")) %>% 
  rownames_to_column(var = "Sample_ID") %>% 
  mutate(Sample_ID = gsub("\\.", "-", Sample_ID)) %>% 
  full_join(., mini_metadata, by="Sample_ID") %>% 
  select(!c(Sample_ID, Sample_Name, Fox_Aspect, Month)) %>% 
  filter(!is.na(Observed)) %>% 
  group_by(DNA_or_cDNA, Snowpack_Layer) %>%
  mutate(Observed = mean(Observed),
         Shannon = mean(Shannon),
         Simpson = mean(Simpson)) %>% 
  distinct_all() %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% 
  row_spec(0, bold=TRUE)
Observed Shannon Simpson DNA_or_cDNA Snowpack_Layer
952 3.2764 0.6659 cDNA TOP
2719 5.7179 0.9837 DNA TOP
453 4.0612 0.9622 cDNA SUP ICE
1739 5.2347 0.9838 DNA SUP ICE

8 - Alpha diversity with rarefaction to the smallest number of reads in any sample

ps.final.rare = rarefy_even_depth(ps.final, 
                                  sample.size = min(sample_sums(ps.final)),
                                  rngseed = TRUE, replace = TRUE, 
                                  trimOTUs = TRUE, verbose = TRUE)
## `set.seed(TRUE)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(TRUE); .Random.seed` for the full vector
## ...
## 25951OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...
plot_richness(ps.final.rare,
              x="Sample_ID", color="Sample_ID",
              measures=c("Observed","Shannon","Simpson")) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Sample_ID") +
  guides(color=guide_legend(ncol=2)) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Month") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

#now without stats, only shapes
plot_richness(ps.final.rare,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

#now without stats, only shapes
plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Fox_Aspect), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))